"""
Multi-Physics Coupling Framework

Implementation of coupled multi-physics solvers that integrate fluid flow,
heat transfer, buoyancy effects, and turbulence modeling. Provides both
segregated and monolithic coupling strategies with advanced convergence
monitoring and stability control.

Key features:
- Segregated coupling (solve each physics sequentially)
- Monolithic coupling (solve all physics simultaneously)
- Under-relaxation for stability
- Convergence monitoring for coupled systems
- Support for various coupling combinations
- Adaptive coupling strategies

Supported Physics Combinations:
1. Flow + Heat Transfer (forced convection)
2. Flow + Heat + Boussinesq (natural convection)
3. Flow + Turbulence (LES turbulent flow)
4. Flow + Heat + Turbulence (turbulent heat transfer)
5. All physics combined (turbulent natural convection)
"""

using LinearAlgebra
using SparseArrays

"""
    MultiPhysicsOptions

Configuration options for multi-physics coupled simulations.

# Fields
- `include_heat_transfer::Bool`: Enable heat transfer coupling
- `include_boussinesq::Bool`: Enable Boussinesq buoyancy effects
- `include_turbulence::Bool`: Enable turbulence modeling
- `coupling_strategy::Symbol`: Coupling approach (:segregated, :monolithic)
- `max_coupling_iterations::Int`: Maximum coupling iterations per time step
- `coupling_tolerance::Float64`: Convergence tolerance for coupling iterations
- `under_relaxation_factors::Dict{Symbol,Float64}`: Relaxation factors for each physics
- `adaptive_relaxation::Bool`: Enable adaptive under-relaxation
- `convergence_monitoring::Bool`: Enable detailed convergence monitoring
- `restart_on_divergence::Bool`: Restart with smaller time step if diverging
"""
Base.@kwdef struct MultiPhysicsOptions
    include_heat_transfer::Bool = true
    include_boussinesq::Bool = true
    include_turbulence::Bool = false
    coupling_strategy::Symbol = :segregated          # :segregated, :monolithic
    max_coupling_iterations::Int = 10
    coupling_tolerance::Float64 = 1e-4
    under_relaxation_factors::Dict{Symbol,Float64} = Dict(
        :momentum => 0.7,
        :pressure => 0.3,
        :temperature => 0.8,
        :turbulence => 0.5
    )
    adaptive_relaxation::Bool = true
    convergence_monitoring::Bool = true
    restart_on_divergence::Bool = true
end

"""
    MultiPhysicsResult

Results from multi-physics coupled simulation.

# Fields
- `u::Array`: x-velocity field
- `v::Array`: y-velocity field
- `w::Union{Array,Nothing}`: z-velocity field (3D only)
- `p::Array`: Pressure field
- `T::Union{Array,Nothing}`: Temperature field (if heat transfer enabled)
- `nu_sgs::Union{Array,Nothing}`: Subgrid-scale viscosity (if turbulence enabled)
- `buoyancy_force::Union{Array,Nothing}`: Buoyancy forces (if Boussinesq enabled)
- `x::Vector{Float64}`: x-coordinates
- `y::Vector{Float64}`: y-coordinates
- `z::Union{Vector{Float64},Nothing}`: z-coordinates (3D only)
- `converged::Bool`: Whether coupled system converged
- `iterations::Int`: Total time steps performed
- `coupling_iterations::Vector{Int}`: Coupling iterations per time step
- `residual_norm::Float64`: Final coupled residual norm
- `solve_time::Float64`: Total solve time (seconds)
- `convergence_history::Vector{Vector{Float64}}`: Residual history for each physics
- `physics_diagnostics::Dict{String,Any}`: Diagnostics from each physics model
- `options::NSOptions`: Base solver options
- `multiphysics_options::MultiPhysicsOptions`: Multi-physics options
- `heat_options::Union{HeatTransferOptions,Nothing}`: Heat transfer options
- `boussinesq_options::Union{BoussinesqOptions,Nothing}`: Boussinesq options
- `turbulence_model::Union{TurbulenceModel,Nothing}`: Turbulence model
"""
Base.@kwdef struct MultiPhysicsResult
    u::Array
    v::Array
    w::Union{Array,Nothing}
    p::Array
    T::Union{Array,Nothing}
    nu_sgs::Union{Array,Nothing}
    buoyancy_force::Union{Array,Nothing}
    x::Vector{Float64}
    y::Vector{Float64}
    z::Union{Vector{Float64},Nothing}
    converged::Bool
    iterations::Int
    coupling_iterations::Vector{Int}
    residual_norm::Float64
    solve_time::Float64
    convergence_history::Vector{Vector{Float64}}
    physics_diagnostics::Dict{String,Any}
    options::NSOptions
    multiphysics_options::MultiPhysicsOptions
    heat_options::Union{HeatTransferOptions,Nothing}
    boussinesq_options::Union{BoussinesqOptions,Nothing}
    turbulence_model::Union{TurbulenceModel,Nothing}
end

"""
    solve_coupled_flow_heat_2d(ns_options::NSOptions, boussinesq_options::BoussinesqOptions, 
                              heat_options::HeatTransferOptions, 
                              multiphysics_options::MultiPhysicsOptions = MultiPhysicsOptions()) -> MultiPhysicsResult

Solve coupled 2D flow and heat transfer with Boussinesq buoyancy effects.

This is the main interface for natural convection problems combining:
- Incompressible Navier-Stokes equations
- Energy equation with convection-diffusion
- Boussinesq approximation for buoyancy

# Arguments
- `ns_options::NSOptions`: Base Navier-Stokes solver options
- `boussinesq_options::BoussinesqOptions`: Boussinesq approximation parameters
- `heat_options::HeatTransferOptions`: Heat transfer parameters
- `multiphysics_options::MultiPhysicsOptions`: Coupling strategy and parameters

# Returns
- `MultiPhysicsResult`: Complete coupled solution with all physics
"""
function solve_coupled_flow_heat_2d(ns_options::NSOptions, boussinesq_options::BoussinesqOptions, 
                                   heat_options::HeatTransferOptions, 
                                   multiphysics_options::MultiPhysicsOptions = MultiPhysicsOptions())
    if ns_options.verbose
        println("🔄 Starting 2D Coupled Flow-Heat Transfer Simulation")
        println("Physics: Flow + Heat Transfer + Boussinesq Buoyancy")
        println("Coupling: $(multiphysics_options.coupling_strategy)")
        println("Ra ≈ $(estimate_rayleigh_number(boussinesq_options))")
    end
    
    start_time = time()
    
    # Setup domain and initialize fields
    x = LinRange(-1, 1, ns_options.n_block * ns_options.N + 1)
    y = LinRange(-1, 1, ns_options.n_block * ns_options.N + 1)
    nx, ny = length(x), length(y)
    
    # Initialize all fields
    u = zeros(nx, ny)
    v = zeros(nx, ny)
    p = zeros(nx, ny)
    T = zeros(nx, ny)
    buoyancy_force = zeros(nx, ny, 2)
    
    # Initialize temperature field
    initialize_temperature_field!(T, boussinesq_options, x, y)
    
    # Apply initial boundary conditions
    apply_heat_boundary_conditions!(T, heat_options.boundary_conditions, x, y)
    
    # Time stepping and coupling
    dt = compute_multiphysics_timestep_2d(ns_options, multiphysics_options, x, y)
    max_steps = min(ns_options.max_steps, 2000)
    
    convergence_history = Vector{Float64}[]
    coupling_iterations = Int[]
    physics_diagnostics = Dict{String,Any}()
    
    # Main coupled time stepping loop
    for step in 1:max_steps
        if multiphysics_options.coupling_strategy == :segregated
            # Segregated coupling approach
            converged_step, coupling_iters, step_residuals = solve_segregated_step_2d!(
                u, v, p, T, buoyancy_force, dt,
                ns_options, boussinesq_options, heat_options, multiphysics_options,
                x, y, step
            )
        else
            # Monolithic coupling approach
            converged_step, coupling_iters, step_residuals = solve_monolithic_step_2d!(
                u, v, p, T, buoyancy_force, dt,
                ns_options, boussinesq_options, heat_options, multiphysics_options,
                x, y, step
            )
        end
        
        push!(coupling_iterations, coupling_iters)
        push!(convergence_history, step_residuals)
        
        # Apply boundary conditions
        apply_coupled_boundary_conditions_2d!(u, v, T, boussinesq_options, heat_options, x, y)
        
        # Convergence check
        overall_residual = compute_coupled_residual_2d(u, v, T)
        
        if ns_options.verbose && step % 100 == 0
            avg_T = mean(T)
            max_buoyancy = maximum(abs.(buoyancy_force[:,:,2]))
            println("Step $step: Overall residual = $(overall_residual), ⟨T⟩ = $(avg_T), max|f_buoy| = $(max_buoyancy)")
        end
        
        if overall_residual < ns_options.tol
            if ns_options.verbose
                println("Coupled system converged after $step steps")
            end
            break
        end
    end
    
    # Compute final diagnostics
    physics_diagnostics["rayleigh_number"] = compute_rayleigh_number(boussinesq_options, maximum(T) - minimum(T), 2.0)
    physics_diagnostics["nusselt_number"] = compute_nusselt_number_2d(T, :left, 2.0, heat_options.thermal_conductivity)
    physics_diagnostics["average_temperature"] = mean(T)
    physics_diagnostics["max_velocity"] = maximum(sqrt.(u.^2 + v.^2))
    physics_diagnostics["total_buoyancy_work"] = sum(abs.(buoyancy_force[:,:,2] .* v))
    
    solve_time = time() - start_time
    final_residual = isempty(convergence_history) ? 0.0 : sum(convergence_history[end])
    
    return MultiPhysicsResult(
        u = u, v = v, w = nothing, p = p, T = T,
        nu_sgs = nothing, buoyancy_force = buoyancy_force,
        x = collect(x), y = collect(y), z = nothing,
        converged = final_residual < ns_options.tol,
        iterations = length(convergence_history),
        coupling_iterations = coupling_iterations,
        residual_norm = final_residual,
        solve_time = solve_time,
        convergence_history = convergence_history,
        physics_diagnostics = physics_diagnostics,
        options = ns_options,
        multiphysics_options = multiphysics_options,
        heat_options = heat_options,
        boussinesq_options = boussinesq_options,
        turbulence_model = nothing
    )
end

"""
    solve_les_turbulent_heat_transfer_2d(ns_options::NSOptions, heat_options::HeatTransferOptions,
                                        turbulence_model::TurbulenceModel,
                                        multiphysics_options::MultiPhysicsOptions = MultiPhysicsOptions()) -> MultiPhysicsResult

Solve coupled 2D LES turbulent flow with heat transfer.

Combines:
- Large Eddy Simulation for turbulent flow
- Energy equation with turbulent mixing
- Enhanced heat transfer due to turbulence

# Arguments
- `ns_options::NSOptions`: Base solver options
- `heat_options::HeatTransferOptions`: Heat transfer parameters  
- `turbulence_model::TurbulenceModel`: LES turbulence model
- `multiphysics_options::MultiPhysicsOptions`: Coupling parameters

# Returns
- `MultiPhysicsResult`: Complete turbulent heat transfer solution
"""
function solve_les_turbulent_heat_transfer_2d(ns_options::NSOptions, heat_options::HeatTransferOptions,
                                             turbulence_model::TurbulenceModel,
                                             multiphysics_options::MultiPhysicsOptions = MultiPhysicsOptions())
    if ns_options.verbose
        println("🌪️🔥 Starting 2D LES Turbulent Heat Transfer Simulation")
        println("Physics: LES Turbulence + Heat Transfer")
        println("Turbulence model: $(typeof(turbulence_model))")
    end
    
    start_time = time()
    
    # Setup domain and initialize fields
    x = LinRange(-1, 1, ns_options.n_block * ns_options.N + 1)
    y = LinRange(-1, 1, ns_options.n_block * ns_options.N + 1)
    nx, ny = length(x), length(y)
    
    # Initialize fields
    u = zeros(nx, ny)
    v = zeros(nx, ny)
    p = zeros(nx, ny)
    T = zeros(nx, ny)
    nu_sgs = zeros(nx, ny)
    strain_rate = zeros(nx, ny)
    
    # Initialize turbulent flow and temperature
    initialize_turbulent_flow_2d!(u, v, x, y, ns_options.nu)
    initialize_temperature_field!(T, nothing, x, y)  # No Boussinesq
    apply_heat_boundary_conditions!(T, heat_options.boundary_conditions, x, y)
    
    # Time stepping
    dt = compute_multiphysics_timestep_2d(ns_options, multiphysics_options, x, y)
    max_steps = min(ns_options.max_steps, 1500)
    
    convergence_history = Vector{Float64}[]
    coupling_iterations = Int[]
    
    # Main turbulent heat transfer loop
    for step in 1:max_steps
        # Solve LES turbulence
        u_filtered = apply_explicit_filter_2d(u, :top_hat)
        v_filtered = apply_explicit_filter_2d(v, :top_hat)
        compute_strain_rate_tensor_2d!(strain_rate, u_filtered, v_filtered, x, y)
        compute_subgrid_viscosity_2d!(nu_sgs, u_filtered, v_filtered, strain_rate, turbulence_model, x, y)
        
        # Solve momentum with turbulent viscosity
        solve_les_momentum_2d!(u, v, p, nu_sgs, dt, ns_options.nu, x, y)
        
        # Solve energy equation with turbulent mixing
        solve_turbulent_heat_equation_2d!(T, u, v, nu_sgs, dt, heat_options, x, y)
        
        # Apply boundary conditions
        apply_les_boundary_conditions_2d!(u, v, ns_options, x, y)
        apply_heat_boundary_conditions!(T, heat_options.boundary_conditions, x, y)
        
        # Convergence monitoring
        residual_flow = compute_les_residual_2d(u, v)
        residual_heat = norm(T) / length(T)
        step_residuals = [residual_flow, residual_heat]
        push!(convergence_history, step_residuals)
        push!(coupling_iterations, 1)  # Direct coupling
        
        overall_residual = sum(step_residuals)
        
        if ns_options.verbose && step % 200 == 0
            avg_nu_sgs = mean(nu_sgs)
            avg_T = mean(T)
            println("Step $step: Flow residual = $(residual_flow), Heat residual = $(residual_heat), ⟨νₛₘ⟩ = $(avg_nu_sgs), ⟨T⟩ = $(avg_T)")
        end
        
        if overall_residual < ns_options.tol
            break
        end
    end
    
    # Compute diagnostics
    turbulence_stats = compute_turbulence_statistics_2d(u, v, nu_sgs, strain_rate)
    physics_diagnostics = merge(turbulence_stats, Dict{String,Any}(
        "average_temperature" => mean(T),
        "nusselt_number" => compute_nusselt_number_2d(T, :left, 2.0, heat_options.thermal_conductivity),
        "turbulent_prandtl_number" => mean(nu_sgs) / (heat_options.thermal_diffusivity + 1e-12)
    ))
    
    solve_time = time() - start_time
    final_residual = isempty(convergence_history) ? 0.0 : sum(convergence_history[end])
    
    return MultiPhysicsResult(
        u = u, v = v, w = nothing, p = p, T = T,
        nu_sgs = nu_sgs, buoyancy_force = nothing,
        x = collect(x), y = collect(y), z = nothing,
        converged = final_residual < ns_options.tol,
        iterations = length(convergence_history),
        coupling_iterations = coupling_iterations,
        residual_norm = final_residual,
        solve_time = solve_time,
        convergence_history = convergence_history,
        physics_diagnostics = physics_diagnostics,
        options = ns_options,
        multiphysics_options = multiphysics_options,
        heat_options = heat_options,
        boussinesq_options = nothing,
        turbulence_model = turbulence_model
    )
end

"""
    solve_full_multiphysics_2d(ns_options::NSOptions, boussinesq_options::BoussinesqOptions,
                              heat_options::HeatTransferOptions, turbulence_model::TurbulenceModel,
                              multiphysics_options::MultiPhysicsOptions = MultiPhysicsOptions()) -> MultiPhysicsResult

Solve full 2D multi-physics problem with all effects combined:
- Turbulent flow (LES)
- Heat transfer with turbulent mixing
- Buoyancy effects (Boussinesq approximation)

This represents the most complete physics combination available.
"""
function solve_full_multiphysics_2d(ns_options::NSOptions, boussinesq_options::BoussinesqOptions,
                                   heat_options::HeatTransferOptions, turbulence_model::TurbulenceModel,
                                   multiphysics_options::MultiPhysicsOptions = MultiPhysicsOptions())
    if ns_options.verbose
        println("🌪️🔥🌡️  Starting Full 2D Multi-Physics Simulation")
        println("Physics: LES Turbulence + Heat Transfer + Boussinesq Buoyancy")
        println("This is the most computationally intensive simulation available")
    end
    
    start_time = time()
    
    # Setup domain
    x = LinRange(-1, 1, ns_options.n_block * ns_options.N + 1)
    y = LinRange(-1, 1, ns_options.n_block * ns_options.N + 1)
    nx, ny = length(x), length(y)
    
    # Initialize all fields
    u = zeros(nx, ny)
    v = zeros(nx, ny)
    p = zeros(nx, ny)
    T = zeros(nx, ny)
    nu_sgs = zeros(nx, ny)
    strain_rate = zeros(nx, ny)
    buoyancy_force = zeros(nx, ny, 2)
    
    # Initialize with turbulent natural convection conditions
    initialize_turbulent_flow_2d!(u, v, x, y, ns_options.nu)
    initialize_temperature_field!(T, boussinesq_options, x, y)
    apply_heat_boundary_conditions!(T, heat_options.boundary_conditions, x, y)
    
    # Time stepping (more conservative for full physics)
    dt = 0.5 * compute_multiphysics_timestep_2d(ns_options, multiphysics_options, x, y)
    max_steps = min(ns_options.max_steps, 1000)  # Reduced for computational cost
    
    convergence_history = Vector{Float64}[]
    coupling_iterations = Int[]
    
    # Full multi-physics time stepping loop
    for step in 1:max_steps
        # Step 1: Solve turbulence
        u_filtered = apply_explicit_filter_2d(u, :top_hat)
        v_filtered = apply_explicit_filter_2d(v, :top_hat)
        compute_strain_rate_tensor_2d!(strain_rate, u_filtered, v_filtered, x, y)
        compute_subgrid_viscosity_2d!(nu_sgs, u_filtered, v_filtered, strain_rate, turbulence_model, x, y)
        
        # Step 2: Compute buoyancy forces
        compute_buoyancy_force!(buoyancy_force, T, 
                               boussinesq_options.density_reference,
                               boussinesq_options.thermal_expansion,
                               boussinesq_options.gravitational_acceleration,
                               boussinesq_options.reference_temperature)
        
        # Step 3: Solve momentum with turbulence and buoyancy
        solve_full_momentum_2d!(u, v, p, nu_sgs, buoyancy_force, dt, ns_options.nu, x, y)
        
        # Step 4: Solve energy equation with turbulent mixing
        solve_turbulent_heat_equation_2d!(T, u, v, nu_sgs, dt, heat_options, x, y)
        
        # Step 5: Apply all boundary conditions
        apply_les_boundary_conditions_2d!(u, v, ns_options, x, y)
        apply_heat_boundary_conditions!(T, heat_options.boundary_conditions, x, y)
        
        # Convergence monitoring
        residual_flow = compute_les_residual_2d(u, v)
        residual_heat = norm(T) / length(T)
        residual_buoyancy = norm(buoyancy_force) / length(buoyancy_force)
        step_residuals = [residual_flow, residual_heat, residual_buoyancy]
        push!(convergence_history, step_residuals)
        push!(coupling_iterations, 1)  # Fully coupled approach
        
        overall_residual = sum(step_residuals)
        
        if ns_options.verbose && step % 100 == 0
            avg_nu_sgs = mean(nu_sgs)
            avg_T = mean(T)
            max_buoyancy = maximum(abs.(buoyancy_force[:,:,2]))
            println("Step $step: Overall residual = $(overall_residual), ⟨νₛₘ⟩ = $(avg_nu_sgs), ⟨T⟩ = $(avg_T), max|f_b| = $(max_buoyancy)")
        end
        
        if overall_residual < ns_options.tol
            break
        end
    end
    
    # Comprehensive diagnostics
    turbulence_stats = compute_turbulence_statistics_2d(u, v, nu_sgs, strain_rate)
    physics_diagnostics = merge(turbulence_stats, Dict{String,Any}(
        "rayleigh_number" => compute_rayleigh_number(boussinesq_options, maximum(T) - minimum(T), 2.0),
        "nusselt_number" => compute_nusselt_number_2d(T, :left, 2.0, heat_options.thermal_conductivity),
        "average_temperature" => mean(T),
        "turbulent_prandtl_number" => mean(nu_sgs) / (heat_options.thermal_diffusivity + 1e-12),
        "buoyancy_reynolds_number" => maximum(sqrt.(u.^2 + v.^2)) * 2.0 / ns_options.nu
    ))
    
    solve_time = time() - start_time
    final_residual = isempty(convergence_history) ? 0.0 : sum(convergence_history[end])
    
    return MultiPhysicsResult(
        u = u, v = v, w = nothing, p = p, T = T,
        nu_sgs = nu_sgs, buoyancy_force = buoyancy_force,
        x = collect(x), y = collect(y), z = nothing,
        converged = final_residual < ns_options.tol,
        iterations = length(convergence_history),
        coupling_iterations = coupling_iterations,
        residual_norm = final_residual,
        solve_time = solve_time,
        convergence_history = convergence_history,
        physics_diagnostics = physics_diagnostics,
        options = ns_options,
        multiphysics_options = multiphysics_options,
        heat_options = heat_options,
        boussinesq_options = boussinesq_options,
        turbulence_model = turbulence_model
    )
end

# Segregated coupling implementation

"""
    solve_segregated_step_2d!(u, v, p, T, f_buoy, dt, ns_opt, bouss_opt, heat_opt, mp_opt, x, y, step)

Solve one time step using segregated coupling approach.
"""
function solve_segregated_step_2d!(u::Array{Float64,2}, v::Array{Float64,2}, p::Array{Float64,2}, 
                                  T::Array{Float64,2}, f_buoy::Array{Float64,3}, dt::Float64,
                                  ns_opt::NSOptions, bouss_opt::BoussinesqOptions, 
                                  heat_opt::HeatTransferOptions, mp_opt::MultiPhysicsOptions,
                                  x::AbstractVector, y::AbstractVector, step::Int)
    
    convergence_residuals = Float64[]
    coupling_iter = 0
    
    # Store old values for convergence checking
    u_old = copy(u)
    v_old = copy(v)
    T_old = copy(T)
    
    # Segregated coupling iterations
    for iter in 1:mp_opt.max_coupling_iterations
        coupling_iter = iter
        
        # Step 1: Update buoyancy forces based on current temperature
        compute_buoyancy_force!(f_buoy, T, 
                               bouss_opt.density_reference,
                               bouss_opt.thermal_expansion,
                               bouss_opt.gravitational_acceleration,
                               bouss_opt.reference_temperature)
        
        # Step 2: Solve momentum equations with buoyancy
        α_mom = mp_opt.under_relaxation_factors[:momentum]
        u_new = copy(u)
        v_new = copy(v)
        
        solve_momentum_with_buoyancy!(u_new, v_new, p, f_buoy, dt, ns_opt.nu, x, y)
        
        # Apply under-relaxation
        u .= α_mom * u_new + (1 - α_mom) * u
        v .= α_mom * v_new + (1 - α_mom) * v
        
        # Step 3: Solve energy equation with updated velocity
        α_temp = mp_opt.under_relaxation_factors[:temperature]
        T_new = copy(T)
        
        thermal_diffusivity = bouss_opt.prandtl_number * ns_opt.nu
        solve_energy_equation!(T_new, T, u, v, dt, thermal_diffusivity, x, y)
        
        # Apply under-relaxation
        T .= α_temp * T_new + (1 - α_temp) * T
        
        # Check coupling convergence
        residual_u = norm(u - u_old) / (norm(u) + 1e-12)
        residual_v = norm(v - v_old) / (norm(v) + 1e-12)
        residual_T = norm(T - T_old) / (norm(T) + 1e-12)
        
        coupling_residual = max(residual_u, residual_v, residual_T)
        push!(convergence_residuals, coupling_residual)
        
        if mp_opt.convergence_monitoring && ns_opt.verbose && step % 100 == 0
            println("  Coupling iter $iter: Residual = $(coupling_residual)")
        end
        
        if coupling_residual < mp_opt.coupling_tolerance
            break
        end
        
        # Update old values for next iteration
        u_old .= u
        v_old .= v  
        T_old .= T
    end
    
    converged = coupling_iter < mp_opt.max_coupling_iterations
    return converged, coupling_iter, convergence_residuals
end

"""
    solve_monolithic_step_2d!(u, v, p, T, f_buoy, dt, ns_opt, bouss_opt, heat_opt, mp_opt, x, y, step)

Solve one time step using monolithic coupling approach (simplified implementation).
"""
function solve_monolithic_step_2d!(u::Array{Float64,2}, v::Array{Float64,2}, p::Array{Float64,2}, 
                                   T::Array{Float64,2}, f_buoy::Array{Float64,3}, dt::Float64,
                                   ns_opt::NSOptions, bouss_opt::BoussinesqOptions, 
                                   heat_opt::HeatTransferOptions, mp_opt::MultiPhysicsOptions,
                                   x::AbstractVector, y::AbstractVector, step::Int)
    
    # For monolithic coupling, we solve all equations simultaneously
    # This is a simplified implementation - full monolithic would require
    # assembling a large coupled system matrix
    
    # Use smaller sub-time steps with explicit coupling
    n_substeps = 4
    dt_sub = dt / n_substeps
    
    convergence_residuals = Float64[]
    
    for substep in 1:n_substeps
        # Compute buoyancy based on current temperature
        compute_buoyancy_force!(f_buoy, T, 
                               bouss_opt.density_reference,
                               bouss_opt.thermal_expansion,
                               bouss_opt.gravitational_acceleration,
                               bouss_opt.reference_temperature)
        
        # Solve momentum and energy simultaneously with smaller time step
        u_old = copy(u)
        v_old = copy(v)
        T_old = copy(T)
        
        solve_momentum_with_buoyancy!(u, v, p, f_buoy, dt_sub, ns_opt.nu, x, y)
        
        thermal_diffusivity = bouss_opt.prandtl_number * ns_opt.nu
        solve_energy_equation!(T, T_old, u, v, dt_sub, thermal_diffusivity, x, y)
        
        # Compute residual for this substep
        residual = max(
            norm(u - u_old) / (norm(u) + 1e-12),
            norm(v - v_old) / (norm(v) + 1e-12),
            norm(T - T_old) / (norm(T) + 1e-12)
        )
        push!(convergence_residuals, residual)
    end
    
    return true, n_substeps, convergence_residuals  # Always "converged" for monolithic
end

# Heat transfer with turbulent mixing

"""
    solve_turbulent_heat_equation_2d!(T_new, T_old, u, v, nu_sgs, dt, heat_options, x, y)

Solve energy equation with turbulent thermal diffusivity.
"""
function solve_turbulent_heat_equation_2d!(T_new::Array{Float64,2}, u::Array{Float64,2}, v::Array{Float64,2},
                                          nu_sgs::Array{Float64,2}, dt::Float64, heat_options::HeatTransferOptions,
                                          x::AbstractVector, y::AbstractVector)
    nx, ny = size(T_new)
    dx = length(x) > 1 ? abs(x[2] - x[1]) : 1.0
    dy = length(y) > 1 ? abs(y[2] - y[1]) : 1.0
    
    T_old = copy(T_new)
    
    # Turbulent Prandtl number (typically 0.7-0.9)
    Pr_t = 0.8
    
    for i in 2:(nx-1), j in 2:(ny-1)
        # Enhanced thermal diffusivity: α_eff = α + νₛₘ/Prₜ
        alpha_molecular = heat_options.thermal_diffusivity
        alpha_turbulent = nu_sgs[i,j] / Pr_t
        alpha_effective = alpha_molecular + alpha_turbulent
        
        # Convection terms
        dTdx = (T_old[i+1,j] - T_old[i-1,j]) / (2*dx)
        dTdy = (T_old[i,j+1] - T_old[i,j-1]) / (2*dy)
        convection = u[i,j] * dTdx + v[i,j] * dTdy
        
        # Enhanced diffusion terms
        d2Tdx2 = (T_old[i+1,j] - 2*T_old[i,j] + T_old[i-1,j]) / dx^2
        d2Tdy2 = (T_old[i,j+1] - 2*T_old[i,j] + T_old[i,j-1]) / dy^2
        diffusion = alpha_effective * (d2Tdx2 + d2Tdy2)
        
        # Heat source
        heat_source = 0.0
        if heat_options.heat_source !== nothing
            heat_source = heat_options.heat_source(x[i], y[j], 0.0, 0.0) / (heat_options.density * heat_options.heat_capacity)
        end
        
        # Update temperature with turbulent mixing
        T_new[i,j] = T_old[i,j] + dt * (-convection + diffusion + heat_source)
    end
end

"""
    solve_full_momentum_2d!(u, v, p, nu_sgs, f_buoy, dt, nu, x, y)

Solve momentum equations with both turbulence and buoyancy effects.
"""
function solve_full_momentum_2d!(u::Array{Float64,2}, v::Array{Float64,2}, p::Array{Float64,2},
                                nu_sgs::Array{Float64,2}, f_buoy::Array{Float64,3},
                                dt::Float64, nu::Float64, x::AbstractVector, y::AbstractVector)
    nx, ny = size(u)
    dx = length(x) > 1 ? abs(x[2] - x[1]) : 1.0
    dy = length(y) > 1 ? abs(y[2] - y[1]) : 1.0
    
    u_old = copy(u)
    v_old = copy(v)
    
    for i in 2:(nx-1), j in 2:(ny-1)
        # Total effective viscosity: molecular + subgrid-scale
        nu_effective = nu + nu_sgs[i,j]
        
        # Convective terms
        dudt_conv = -u_old[i,j] * (u_old[i+1,j] - u_old[i-1,j])/(2*dx) - 
                    v_old[i,j] * (u_old[i,j+1] - u_old[i,j-1])/(2*dy)
                    
        dvdt_conv = -u_old[i,j] * (v_old[i+1,j] - v_old[i-1,j])/(2*dx) - 
                    v_old[i,j] * (v_old[i,j+1] - v_old[i,j-1])/(2*dy)
        
        # Pressure gradients
        dpdx = (p[i+1,j] - p[i-1,j])/(2*dx)
        dpdy = (p[i,j+1] - p[i,j-1])/(2*dy)
        
        # Enhanced viscous terms
        d2udx2 = (u_old[i+1,j] - 2*u_old[i,j] + u_old[i-1,j])/(dx^2)
        d2udy2 = (u_old[i,j+1] - 2*u_old[i,j] + u_old[i,j-1])/(dy^2)
        d2vdx2 = (v_old[i+1,j] - 2*v_old[i,j] + v_old[i-1,j])/(dx^2)
        d2vdy2 = (v_old[i,j+1] - 2*v_old[i,j] + v_old[i,j-1])/(dy^2)
        
        # Update with turbulence and buoyancy
        u[i,j] = u_old[i,j] + dt * (dudt_conv - dpdx + nu_effective*(d2udx2 + d2udy2) + f_buoy[i,j,1])
        v[i,j] = v_old[i,j] + dt * (dvdt_conv - dpdy + nu_effective*(d2vdx2 + d2vdy2) + f_buoy[i,j,2])
    end
end

# Utility and helper functions

"""
    initialize_temperature_field!(T, boussinesq_options, x, y)

Initialize temperature field for multi-physics simulations.
"""
function initialize_temperature_field!(T::Array{Float64,2}, boussinesq_options::Union{BoussinesqOptions,Nothing},
                                      x::AbstractVector, y::AbstractVector)
    nx, ny = size(T)
    
    if boussinesq_options !== nothing
        # Temperature field for natural convection
        T_ref = boussinesq_options.reference_temperature
        T_hot = T_ref + 10.0
        T_cold = T_ref - 10.0
        
        for (i, xi) in enumerate(x), (j, yj) in enumerate(y)
            # Linear temperature profile
            T[i,j] = T_cold + (T_hot - T_cold) * (xi + 1) / 2
        end
    else
        # Temperature field for forced convection
        for (i, xi) in enumerate(x), (j, yj) in enumerate(y)
            if i == 1  # Left wall heated
                T[i,j] = 350.0  # Hot wall
            elseif i == nx  # Right wall cooled
                T[i,j] = 300.0  # Cold wall
            else
                T[i,j] = 325.0  # Initial bulk temperature
            end
        end
    end
end

"""
    apply_coupled_boundary_conditions_2d!(u, v, T, bouss_opt, heat_opt, x, y)

Apply boundary conditions for coupled flow-heat system.
"""
function apply_coupled_boundary_conditions_2d!(u::Array{Float64,2}, v::Array{Float64,2}, T::Array{Float64,2},
                                              bouss_opt::BoussinesqOptions, heat_opt::HeatTransferOptions,
                                              x::AbstractVector, y::AbstractVector)
    # Flow boundary conditions (no-slip walls)
    apply_boussinesq_boundary_conditions!(u, v, T, bouss_opt, x, y)
    
    # Additional heat transfer boundary conditions
    apply_heat_boundary_conditions!(T, heat_opt.boundary_conditions, x, y)
end

"""
    compute_coupled_residual_2d(u, v, T) -> Float64

Compute overall residual for coupled system.
"""
function compute_coupled_residual_2d(u::Array{Float64,2}, v::Array{Float64,2}, T::Array{Float64,2})
    residual_flow = sqrt(sum(u.^2) + sum(v.^2)) / (length(u) + length(v))
    residual_heat = sqrt(sum(T.^2)) / length(T)
    
    # Weighted combination
    return 0.5 * residual_flow + 0.5 * residual_heat
end

"""
    compute_multiphysics_timestep_2d(ns_options, mp_options, x, y) -> Float64

Compute stable time step for multi-physics simulation.
"""
function compute_multiphysics_timestep_2d(ns_options::NSOptions, mp_options::MultiPhysicsOptions,
                                         x::AbstractVector, y::AbstractVector)
    dx = length(x) > 1 ? abs(x[2] - x[1]) : 1.0
    dy = length(y) > 1 ? abs(y[2] - y[1]) : 1.0
    h_min = min(dx, dy)
    
    # Base time step from flow solver
    dt_flow = ns_options.cfl * h_min^2 / (2 * ns_options.nu)
    
    # Additional restrictions for multi-physics
    dt_restrictions = [dt_flow]
    
    if mp_options.include_heat_transfer
        # Heat transfer stability
        alpha_typical = 2.2e-5  # Typical thermal diffusivity
        dt_heat = 0.4 * h_min^2 / (2 * alpha_typical)
        push!(dt_restrictions, dt_heat)
    end
    
    if mp_options.include_turbulence
        # More conservative time step for turbulence
        dt_turbulent = 0.3 * dt_flow
        push!(dt_restrictions, dt_turbulent)
    end
    
    if mp_options.coupling_strategy == :segregated
        # Segregated coupling can be less stable
        dt_coupling = 0.5 * minimum(dt_restrictions)
        push!(dt_restrictions, dt_coupling)
    end
    
    return minimum(dt_restrictions)
end

"""
    estimate_rayleigh_number(boussinesq_options) -> Float64

Estimate Rayleigh number from Boussinesq parameters.
"""
function estimate_rayleigh_number(boussinesq_options::BoussinesqOptions)
    g_mag = norm(boussinesq_options.gravitational_acceleration)
    β = boussinesq_options.thermal_expansion
    ΔT = 20.0  # Typical temperature difference
    L = 2.0    # Characteristic length
    
    # Assume typical values for missing properties
    ν = 1.5e-5  # Kinematic viscosity (air)
    α = 2.2e-5  # Thermal diffusivity (air)
    
    Ra = g_mag * β * ΔT * L^3 / (ν * α)
    return Ra
end

"""
    monitor_coupling_convergence(residuals, tolerances) -> Bool

Monitor convergence of coupling iterations with detailed diagnostics.
"""
function monitor_coupling_convergence(residuals::Vector{Vector{Float64}}, tolerances::Vector{Float64})
    if isempty(residuals)
        return false
    end
    
    current_residuals = residuals[end]
    n_physics = length(current_residuals)
    
    if length(tolerances) != n_physics
        @warn "Tolerance vector length mismatch"
        return false
    end
    
    # Check if each physics has converged
    converged_physics = [current_residuals[i] < tolerances[i] for i in 1:n_physics]
    
    # All physics must converge
    return all(converged_physics)
end